Zero Energy Modes and Gate- Tunable Gap in Graphene on hexagonal Boron Nitride 
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In this Letter, we derive an effective theory of graphene on a hexagonal Boron Nitride (h-BN) 
substrate. We show that the h-BN substrate generically opens a spectral gap in graphene despite 
the lattice mismatch. The origin of that gap is particularly intuitive in the regime of strong coupling 
between graphene and its substrate, when the low-energy physics is determined by the topology of a 
network of zero energy modes. For twisted graphene bilayers, where inversion symmetry is present, 
this network percolates through the system and the spectrum is gap less. The breaking of that 
symmetry by h-BN causes the zero energy modes to close into rings. The eigenstates of these rings 
hybridize into flat bands with gaps in between. The size of this band gap can be tuned by a gate 
voltage and it can reach the order of magnitude needed to confine electrons at room temperature. 

PACS numbers: 71.27. +a,73.20.Hb,75.30.Hx 



Graphene is a two-dimensional semimetal with low- 
energy excitations that obey the massless Dirac equa- 
tion [1-3]. As most applications in electronics require a 
bandgap, much effort has been exerted to find ways of in- 
ducing a gap in the electronic spectrum of graphene. One 
possible route is to use hexagonal boron nitride (h-BN) 
substrates 0, 0|, which lack sublattice inversion sym- 
metry. If inherited by the graphene layer, this broken 
symmetry leads to the opening of a gap in the spec- 
trum, which is described by a mass in the Dirac model. 
First principles calculations for graphene supported by 
a perfectly lattice matched h-BN substrate predicted a 
bandgap on the order 50 meV jj|. Experiments, however, 
have not observed any clear indication of a gap 0, [?J . 

Subsequent theory 0, [gj identified the reason for this 
discrepancy: the lattice constants of h-BN and graphene 
differ by about 1.8% and for multi-layered h-BN sub- 
strates it is energetically unfavorable for graphene and 
its substrate to conform their lattice constants [8|- The 
result is Moire patterns of varying local lattice alignment, 
as illustrated in Fig.[IJ which have been observed in scan- 
ning tunneling microscopy [10]. Evidently, the h-BN sub- 
strate breaks sublattice symmetry differently in different 
regions of the Moire pattern. Density functional theory 
(DFT) calculations have shown a tendency for the sub- 
lattice symmetry breaking to be compensated between 
different regions of the Moire pattern such that the sym- 
metry is almost restored after spatial average [8[. This 
has motivated the proposal of effective Dirac models of 
graphene on h-BN with a mass term that has vanish- 
ing integral, such that sublattice inversion symmetry is 
restored on (spatial) average. On the basis of those mod- 
els, it was argued that a band gap in graphene on h-BN is 
absent consistently with the existing experimental 

data 0,0]. 

In this Letter, we derive an effective theory for 
graphene on h-BN based on a bilayer model that has been 
successfully applied to twisted graphene bilayers [Tll-[l3| . 



Our theory is formulated for a single-layer of graphene 
and it accounts for the coupling to the substrate by a 
mass term and effective potentials that oscillate with the 
period of the Moire pattern. We find that graphene sup- 
ported by a hexagonal substrate generically develops a 
gap in the spectrum. In particular, the emergence of a 
spectral gap is not precluded by an average sublattice 
symmetry. A gap is avoided only by additional sym- 
metries, such as in twisted graphene bilayers, where the 
Dirac points are topologically protected by a combina- 
tion of space inversion and time reversal symmetry [lij . 
In the regime of strong coupling between graphene and 
its substrate a particularly intuitive picture of that gap 
formation emerges: the oscillatory mass in our effective 
theory then defines one-dimensional modes that are topo- 
logically protected and gapless for a large Moire period. 
In the presence of space inversion symmetry, such as in 
twisted graphene bilayers, these modes form a network 
that percolates through the system, corresponding to a 
metallic state. When space inversion symmetry is bro- 
ken, for instance by h-BN, this network breaks up into 
isolated rings of ID-modes. The states in these rings 
hybridize exponentially weakly and form narrow bands 
with large gaps set by the level spacing in the rings. 
Finally, we show that the spectral gap of graphene on 




FIG. 1: Two layer system made of graphene (top layer) and 
h-BN (bottom layer) with a lattice mismatch (exaggerated in 
the figure). Red line: Moire unit cell of the system. 
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h-BN can be controlled by the application of a perpen- 
dicular electric field. Using parameters fitted to experi- 
ments and DFT calculations, we find that the gap can be 
tuned to reach the order of magnitude needed to confine 
charge carriers at room temperature, a key requirement 
for electronics applications. 

Model: We base our analysis of graphene with an h- 
BN substrate on a tight-binding model of two coupled 
honeycomb lattices with parameters fitted to experiment 
[l5| . The two layers have a lattice mismatch of S « 1.8% 
(cf. Fig. [T]) and we allow in addition a rotational mis- 
alignment by angle 0. In the two- layer basis of graphene 
and h-BN, the electronic Hamiltonian is 

In the limit <C 1 and S <C 1, where a long- wavelength 
description is appropriate, the isolated graphene layer 
can be described by a Dirac model Hamiltonian H g = 
vp • 5], and the h-BN layer is similarly described by 
Hbn = ^bnP • 5] + mv 2 a z + V, where v and ^bn are 
the carrier velocities in graphene and h-BN, respectively 
and p is the momentum relative to the Dirac points. 
H g and Hbn each act on 4 component spinors of the 
form (^a,+, ^b,+, ^a,-, ^b,-), with A, B labeling the 
two different sublattices of the honeycomb lattice and 
± denoting the two different band structure valleys. 
5] = (t z (j Xi <j y ) is a vector of Pauli matrices acting on the 
A/B sublattice basis (through the pseudospin a) and the 
valley spin (r). The mass m and the interlayer bias V 
account for the different on-site potentials of the Boron 
(the A sublattice of h-BN) and the Nitrogen (B-)atoms. 
DFT calculations indicate mv 2 « 2.3 eV and V w 0.8 eV 
Q . The bias V can be tuned by a perpendicular electric 
field. 

Following Ref. [13] . we take the interlayer coupling of 
Eq. (1) at long- wavelengths to be of the form [111. Il3| 

Here, H{ nt is written explicitly in the A, B sublattice 
basis, 7 « 0.3 eV is the hopping energy to the sub- 
strate, ( parametrizes a sublattice asymmetry due to 
structural differences between different regions of the 
Moire pattern [T^j, and </> n = 27rn/3 is a phase that de- 
pends on the index n = 0,1,2, which labels the three 
corners of the graphene Brillouin zone corresponding to 
a given valley. Those points have wavevectors K n = 
Rfyn + 7r/6)a47r/3v // 3a, where a is a unit vector along an 
A-B bond, a the bond length, and R(cp) is a rotation by 
angle ip. The wavevectors 5K n are the differences of K n 
and their counterparts in the closest valley of the h-BN 
layer. They are shorter than K n by a factor SK n /K n = 
v / £ 2 -2(l + £)(cos0-l)/(l + S) and rotated with re- 
spect to K n by angle = arctan [sin 0/(1 + S — cos0)]. 



We neglect commensuration effects, which are small for 
5,0 < 1 M- 

Effective single-layer theory: Integrating out the elec- 
trons in the h-BN layer, we arrive at an effective Hamil- 
tonian H* s (lj) = H g + 5Hg S (uj) for the graphene layer 
with (we set h = 1) 

SHf(uj) = H int (uj - Hb^H^. (3) 

The mass term m dominates the Hamiltonian Hbn for 
all wavevectors where the employed Dirac model holds. 
At those momenta, p <C K ni we may set ^bn = to 
a good approximation, resul ting in an effective Hamilto- 
nian which is local in space [20| , 

SHf^ = — Hi nt ( n ? ) ^ int ' ( 4 ) 

s uj — V + mv z \ 1 J 

where n parametrizes the inversion symmetry breaking 
through the h-BN substrate, 

r] = (u - V + mv 2 )/(uj - V - mv 2 ). (5) 

The effective Hamiltonian (j4]) can be parametrized in 
terms of effective potentials that oscillate in space with 
the periodicity of the bilayer Moire pattern, 

5Hf = V eff (r) + vveY* ■ A eff (r) + m e ^(r)v 2 a z . (6) 

For perfect rotational alignment, = $ = 0, the effec- 
tive vector potential A eff may be gauged away. In the 
more general case <f> ^ the vector potential generates a 
pseudo-magnetic field and it satisfies the Coulomb gauge 
condition V- A eff = at $ = it/ 2. The mass term m eff (r) 
breaks the sublattice exchange symmetry locally and it 
opens a local gap in the spectrum wherever it exceeds 
1/vL, such that the wavefunctions are localized on the 
length scale L of the Moire pattern. 

A global gap in the spectrum is nevertheless precluded 
when the effective theory Eq. (|4|) is invariant under P = 
o- x t x R(-k) [2l|, that is sublattice exchange a x coupled 
with point reflection R(tt) and valley exchange t x , such 
that the underlying lattice model has inversion and time 
reversal symmetry, as for twisted graphene bilayer s. In 
the following, we analyze the spectral gap of graphene on 
h-BN in the absence of that symmetry, which is explicitly 
broken by the inequivalence of the B and N sites. 

Perturbation theory: We start our analysis of the effec- 
tive theory Eq. (j4]) by a perturbative calculation, valid 
for weak coupling j 2 /\cj — V ±mv 2 \ <C vSK. For a lattice 
mismatch of 5 ~ 1.8% one has v5K 0.22 eV, where the 
lower bound corresponds to = 0. Since 7 w 0.3 eV and 
in the absence of an external bias \V =b m\ > 1.5 eV, the 
system is well in this perturbative regime at low energies 
(M < \V±m\ eV). 

When 77 = 1 the model has inversion symmetry P in 
addition to^time reversal invariance and topological ar- 
guments [14] require the presence of at least two Dirac 
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FIG. 2: Energy spectra in the Moire Brillouin zone for strong 



coupling e — L~f 2 /[(V 



^> 1. (Left) rj = 1, as in 



twisted graphene bilayers (corresponding to the triangular 
network of zero energy modes shown in Fig. [3] a), with a 
metallic spectrum. (Right) (" = 1 and 77 = —0.5 (correspond- 
ing to isolated rings of low-energy modes, as shown in Fig. [3] 
b-d), with narrow bands separated by large gaps (red circle). 
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points (or arcs) [22| in the Moire Brillouin zone. This is 
the situation in twisted graphene bilayers. 

When the sublattice symmetry of the substrate is bro- 
ken (77 ^ 1) such as through h-BN, a gap is not precluded 
by symmetry anymore. In the case ( = 1 it turns out 
that the spatial average (m eff ) = j drm eff (r) vanishes 
and correspondingly no gap is found to leading order in 
the effective potentials. Unlike previously assumed [sl@|, 
this restoration of a symmetry on (spatial) average, how- 
ever, does not suppress the gap in the spectrum entirely. 
A gap does appear at third order in 5Hf [23|: 



77(1 - 77)(2cos2$ - 1) 



r 



Slv 2 5K 2 (V -mv 2 ) 2 



• (7) 



In the more general situation, when £ 7^ 1, the spatial 
average of the mass term (m eff ) is non-zero, and the band 
gap appears already at leading order in perturbation, 



(1 -C 2 ) 



7 



3(V — mv 2 ) 



(8) 



A recent DFT calculation [8] predicted A « 4meV in 
the absence of an external perpendicular field, when 77 & 
—0.5. Assuming 7 = 0.3 eV, we estimate |1 — £ 2 | ~ 0.14 
from Eq. (j8]). The relative magnitude of the local gaps in 
various Moire regions found in the DFT calculation [U 
indicate ( > 1, so we conclude that £ « 1.07. A direct 
fit to the local gaps of Ref. [H in the AA-, AB-, and BA- 
stacked regions, respectively, yields ( = 1.19, r] = — 0.72, 
and 7 = 0.25 eV. One possible reason for the discrepancy 
is higher Fourier harmonics of 5H e ^ that we neglect. 

Eq. (j8]) predicts that the spectral gap of graphene on 
h-BN can be substantially increased by application of a 
perpendicular electric field that decreases the B/N on- 
site energies V =b mv 2 . As the gap increases, perturba- 
tion theory eventually breaks down, and the gap has a 
crossover to a nonperturbative regime. 

Nonperturbative regime: In the nonperturbative limit 
7 2 /|o; — V ±mv 2 \ ^> v5K the spectrum is gapped locally 
in regions where m eff (r) 7^ 0. The low-energy physics of 



FIG. 3: Spatial dependence of the mass term resulting from 
the local lattice misalignment in a bilayer as shown in Fig. 
[TJ The red line indicates the Moire unit cell, and the zb signs 
specify the sign of the mass, a) 77 = 1: metallic state, with 
percolating zeros of ra eff ; b) 77 = 0.9; c) 77 = 0.5; d) 77 = 
—0.5: the network of zeros breaks up into rings. In the strong 
coupling regime, those rings contain fully localized states. 



the system is then dominated by one-dimensional modes 
along the zeros of m eff 17|, [l8[. In the absence of in- 
tervalley scattering and in the limit of a large Moire 
size L — ^ 00 these modes are guaranteed to be gapless 
by topological arguments [19j : the lines with m eff = 
separate regions with effective masses of opposite signs, 
as indicated by the "+" and "-" signs in Fig. [3] The 
massive, single- valley Dirac Hamiltonian H = g(k) • cr, 
where g = v(k x ,k y ,m e ^v) and a = (a x ,a y ,a z ), has a 
topological charge N 3 = J dk x dk y g-(d kx g x d ky g) /47r = 



m eff /(2|m et± |) associated with it (here, g = g/|g|). The 
difference of the charges and A 7 ^ to both sides of a 
line m eff = enforces \N^ — N%\ = 1 zero modes per 
valley [19[. For rotational alignment $ = 0, when A eff 
is a pure gauge and for strong screening of the scalar 
potential F eff , the qualitative low-energy physics thus is 
determined entirely by the topology of the zeros of m eff . 

When 77 = 1, such as in twisted graphene bilayers, the 
zeros of m eff (r) form a triangular network that percolates 
through the entire system, as seen in Fig. [3^. The ex- 
pected pairs of zero modes along these lines provide an 
intuitive explanation for the metallic state required by 
the topological arguments quoted earlier [lij . 

For graphene on h-BN with 77 ^ 1 on the other hand, 
this triangular network of zero energy modes breaks 
up into isolated rings (cf. Fig. [3] b-d). The low-energy 
modes confined to the rings of m eff (r) = form discrete 
states with an energy separation set by the circumference 
C — L of those rings. The states hosted by neighbor- 
ing rings have an overlap that is exponentially small in 
L7 2 / (V — mv 2 )v and form correspondingly narrow (val- 
ley) pairs of bands with gaps A w 2ttv/C in between. Ev- 
idently, in this strong coupling regime, a gap of that or- 
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FIG. 4: Left panels: numerical scaling of the gap A oc 
77, £) as a function of e — L^ 2 /[(V — mv 2 )v] for real- 
istic parameter values at — 0. In weak coupling (s <C 1), 
the predicted crossover from linear scaling (D oc e) at Q 7^ 1 
(a,b), [Eq. ©], to cubic (L> oc s 3 ) at C = 1 (c) [Eq. ©], 
is confirmed. Right: Predicted gap as a function of the in- 
terlayer bias V (tunable by a perpendicular electric field) for 
parameter values from two fits to DFT: 7 = 0.3 eV, n = —0.5, 
C = 1.07 (d) and 7 = 0.25 eV, n = -0.72, C = 1.19 (e) (see 
text). 

der of magnitude will appear regardless of whether (m eff ) 
vanishes or not. This gives a physical reason why an av- 
erage sublattice symmetry cannot prevent the opening of 
a gap in the spectrum. 

For graphene on h-BN the above nonperturbative con- 
siderations do not strictly apply: there is no V such that 
the limit 7 2 /|cj — V =L mv 2 \ ^> vSK is reached for all 
energies uj inside the predicted gap of order 2ttv/C ~ 
vSK > 220 meV. We thus next perform numerical cal- 
culations that use a tight-binding model. Computations 
for the true size of the unit cell at = are challenging. 
We therefore exploit the scale invariance of our theory, 
expressing 

A = D(~i 2 /(V - mv 2 )v5K, 77, () x vSK (9) 

in terms of a function D(e,r]X) that may be evaluated 
for smaller unit cells |24|. The scaling parameter e = 
7 2 l\(y ~ mv 2 )v5K] separates the weak (e <C 1) from the 
strong (e ^> 1) coupling regime. 

In Fig. [3^-c, we plot the scaling function D found from 
tight-binding calculations on a unit cell containing 512 
atoms for typical parameters. When £ = 1, the scaling 
of the gap is cubic in e for weak coupling, crossing over 
to linear (A oc e) behavior for £ 7^ 1, in agreement with 
the perturbative analysis of Eqs. (O and (j8j). Fig. [3]i-e 
shows the gap as a function of V for parameters taken 
from the two above fits to DFT data [25j. Despite the 
uncertainty of the parameters entering our model these 
calculations clearly suggest that it is possible to induce 
gaps in graphene on h-BN on the order of room temper- 
ature. 

Conclusions: We have derived a low-energy theory for 
graphene on hexagonal substrates. Our theory demon- 



strates that a h-BN substrate opens a gap in the spec- 
trum of graphene through a breaking of inversion sym- 
metry even when the sublattice symmetry of graphene is 
restored on (spatial) average. We moreover have shown 
that perpendicular electric fields may be used to enhance 
the predicted gaps up to the scale of room temperature. 
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